home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / VideoToolbox 96.06.15 / VideoToolboxSources / Binomial.c < prev    next >
C/C++ Source or Header  |  1995-07-19  |  8KB  |  276 lines
  1. /*
  2. Binomial.c
  3. Copyright (c) 1990,1991,1992,1994 Denis G. Pelli
  4.  
  5. Various routines that deal with Binomial statistics, including generating random
  6. samples and computing confidence intervals.
  7.  
  8. This file contains two generations of routines that generate confidence
  9. intervals. The original routines BinomialUpperBound and BinomialLowerBound are
  10. old, but are more robust. The new routines are based on a Numerical Recipes in C
  11. routine to compute the incomplete beta function. The new routines add little to
  12. the old routines, and require the Numerical Recipes in C, so I commented them
  13. out. 
  14.  
  15. Also see: ChiSquare.c, Exponential.c, Normal.c, Uniform.c
  16.  
  17. HISTORY:
  18. 1/5/91 Added new binomial routines. These routines may be 
  19.         commented out by setting NEW_BINOMIAL to zero.
  20. 1/7/91 Rewrote BinomialUpperBound & BinomialLowerBound to allow specification of the
  21.         desired confidence of the interval.
  22. 4/24/92    dgp    Added BinomialSampleQuickly().
  23. 4/27/92    dgp    Oops. I just noticed that this file was always including nr.h, which you'd
  24.             only have if you own the Numerical Recipes. I moved the include statement
  25.             down inside the conditional below, where it belongs. 
  26. */
  27. #include "VideoToolbox.h"
  28.  
  29. #define NEW_BINOMIAL 0
  30.  
  31. long BinomialSample(double p,long n)
  32. /*
  33. Returns a a sample from a binomial distribution: number of heads in n flips of
  34. a weighted coin with a probability p of a head on each flip. This is fine if n is small. However,
  35. if n is large, e.g. n>100, then this routine will be rather slow, and you may prefer
  36. to use the more elaborate Numerical Recipes bnldev() routine.
  37. */
  38. {
  39.     long k,i;
  40.  
  41.     k=0;
  42.     for(i=0;i<n;i++) if(p>UniformSample())k++;
  43.     return k;
  44. }
  45.  
  46. int BinomialSampleQuickly(int n)
  47. /*
  48. Does n flips of an unweighted coin and returns the number of heads. Very fast if n<1000. 
  49. For larger n it would be faster to use bnldev().
  50. */
  51. {
  52.     register int i,k=0;
  53.     register short r;
  54.     
  55.     for(i=n;i>=8;i-=8){
  56.         r=randU();    /* Use only the upper 8 bits, which are reputed to be more random. */
  57.         if(r<0)k++;
  58.         r<<=1;
  59.         if(r<0)k++;
  60.         r<<=1;
  61.         if(r<0)k++;
  62.         r<<=1;
  63.         if(r<0)k++;
  64.         r<<=1;
  65.         if(r<0)k++;
  66.         r<<=1;
  67.         if(r<0)k++;
  68.         r<<=1;
  69.         if(r<0)k++;
  70.         r<<=1;
  71.         if(r<0)k++;
  72.         r<<=1;
  73.     }
  74.     if(i>0){
  75.         r=randU();
  76.         for(;i>0;i--){
  77.             if(r<0)k++;
  78.             r<<=1;
  79.         }
  80.     }
  81.     return k;
  82. }
  83.  
  84. /*
  85. BinomialUpperBound and BinomialLowerBound return a P confidence interval for the
  86. underlying binomial probability assumed to have generated the data.
  87.  
  88. The formula is based on a Gaussian approximation, solving for the p that will
  89. put the observed result the right number of standard deviations away from the
  90. mean, plus or minus 0.5, as a "continuity correction".
  91.  
  92. This is the best formula I could find in looking through several statistics
  93. books. However, the result has to be taken with a grain of salt because it is
  94. not possible to produce a binomial confidence interval that will satisfy the
  95. strict definition of a confidence interval, namely one that will have the
  96. specified probability P of containing the unknown but fixed parameter p. That's
  97. because, unlike the Normal distribution, the Binomial distribution is not
  98. translation invariant. 
  99.  
  100. I gave some thought to taking a Bayesian approach, assuming a uniform prior pdf
  101. for p and then computing an a posteriori confidence interval. This can be done,
  102. but it's hard. I concluded that it's pointless because the uniform prior pdf
  103. assumption is usually unwarranted.
  104.  
  105. Even though this confidence interval is not wholly satisfactory (i.e. fundamentally
  106. false) it is useful in practice since it typically behaves similarly to the
  107. Normal case, which is theoretically sound.
  108. */
  109. double BinomialLowerBound(double P,long k,long n)
  110. /*
  111. Arguments are the confidence P, and the number k of heads in n flips of a coin.
  112. The returned value p is the lower end of a 2?-sided P confidence interval on the
  113. underlying probability of a head on a single trial.
  114. */
  115. {
  116.     double right,s,ss,p;
  117.     
  118.     if(k>0 && k<n){
  119.         right=k-0.5;
  120.         s=InverseNormal(sqrt(P));
  121.         ss=s*s;
  122.         p=(right+0.5*ss-s*sqrt(right*(1.0-right/n)+0.25*ss))/(n+ss);
  123.         return p;
  124.     }
  125.     if(k==0) return 0.0;
  126.     if(k==n) return pow(1.0-P,1.0/n);
  127.     return sqrt(-1.0);    /* domain error */
  128. }
  129.  
  130. double BinomialUpperBound(double P,long k,long n)
  131. /*
  132. Arguments are the confidence P, and the number k of heads in n flips of a coin.
  133. The returned value p is the upper end of a 2?-sided P confidence interval on the
  134. underlying probability of a head on a single trial.
  135. */
  136. {
  137.     double right,s,ss,p;
  138.     
  139.     if(k>0 && k<n){
  140.         right=k+0.5;
  141.         s=InverseNormal(sqrt(P));
  142.         ss=s*s;
  143.         p=(right+0.5*ss+s*sqrt(right*(1.0-right/n)+0.25*ss))/(n+ss);
  144.         return p;
  145.     }
  146.     if(k==0) return 1.0-pow(1.0-P,1.0/n);
  147.     if(k==n) return 1.0;
  148.     return sqrt(-1.0);    /* domain error */
  149. }
  150.  
  151. #if NEW_BINOMIAL    /* new routines that require Numerical Recipes in C */
  152.     #include "nr.h"                /* prototype for betai() */
  153.     #if 0
  154.         void main(void)
  155.         /* a quick and dirty driver to test some of these routines */
  156.         {
  157.             double p,pUpper,pLower,P,PUpper,PLower;
  158.             int i,n,k;
  159.             
  160.             Require(0);
  161.             n=10;
  162.             p=.2;
  163.             pLower=InverseBinomial(0.05,21,50);
  164.             printf("%f, Abramowitz & Stegun p. 960 Example 18 say this should be 0.3003.\n",pLower);
  165.             P=0.95;
  166.             for(i=0;i<10;i++){
  167.                 k=BinomialSample( p, n);
  168.                 pLower=BinomialLowerBound(P,k,n);
  169.                 pUpper=BinomialUpperBound(P,k,n);
  170.                 PLower=1.0-Binomial(pLower,k,n);
  171.                 PUpper=Binomial(pUpper,k+1,n);
  172.                 printf("n %4d p %6.4f %6.4f≤%6.4f≤%6.4f  %6.4f< <%6.4f\n"
  173.                     ,n,p,pLower,k/(double)n,pUpper,PLower,PUpper);
  174.                 pLower=InverseBinomial(1.0-sqrt(P),k,n);
  175.                 pUpper=InverseBinomial(sqrt(P),k+1,n);
  176.                 PLower=1.0-Binomial(pLower,k,n);
  177.                 PUpper=Binomial(pUpper,k+1,n);
  178.                 printf("  %4d p %6.4f %6.4f≤%6.4f≤%6.4f  %6.4f< <%6.4f\n"
  179.                     ,n,p,pLower,k/(double)n,pUpper,PLower,PUpper);
  180.             }
  181.         }
  182.     #endif
  183.     
  184.     double Binomial(double p,long k,long n)
  185.     /*
  186.     Returns the probability of k or more heads in n flips, where probability of
  187.     each head is p. This identity appears in Numerical Recipes in C, page 182, and
  188.     in Abramowitz and Stegun, p. 945. Eq.26.5.24.
  189.     */
  190.     {
  191.         if(k>n)return 0.0;
  192.         if(k<=0L)return 1.0;
  193.         if(p==0.0){
  194.             if(k==0)return 1.0;
  195.             else return 0.0;
  196.         }
  197.         if(p==1.0)return 1.0;
  198.         return IncompleteBeta(p,k,n-k+1);
  199.     }
  200.     
  201.     double BinomialPdf(double p,long k,long n)
  202.     /*
  203.     Returns the probability of exactly k heads in n flips, where probability of
  204.     each head is p. 
  205.     I'm not really sure whether this is more computationally efficient than computing
  206.     it directly, from the definition of the binomial distribution, but it is very easy
  207.     to write, and avoids the difficulty of computing the binomial coefficient
  208.     without overflow.
  209.     */
  210.     {
  211.         if(k<0L || k>n)return 0.0;
  212.         if(k==n) return pow(p,n);
  213.         return Binomial(p,k,n)-Binomial(p,k+1,n);
  214.     }
  215.     
  216.     double InverseBinomial(double P,long k,long n)
  217.     /* Returns the P-th quantile for the probability p of a heads,
  218.     given k heads in n flips. 
  219.     pUpper=InverseBinomial(0.975,k+1,n) is a 97.5% confidence upper bound on p
  220.     pLower=InverseBinomial(0.025,k,n) is a 97.5% confidence lower bound on p
  221.     Taken together these bounds form a 97.5%*97.5%=95% confidence interval: [pLower,pUpper].
  222.     */
  223.     {
  224.     #if 0
  225.         return InverseIncompleteBeta(P,k,n-k+1);
  226.     #else
  227.         /*
  228.         This simple-minded bisection routine is slow, but its answer is accurate
  229.         to within ±1e-10. 
  230.         */
  231.         double low=0.0,high=1.0,mid;
  232.         double f;
  233.         int i;
  234.         
  235.         for(i=0;i<30;i++){
  236.             mid=(low+high)/2.0;
  237.             f=Binomial(mid,k,n);
  238.             if(f>P)high=mid;
  239.             else low=mid;
  240.         }
  241.         return (low+high)/2.0;
  242.     #endif
  243.     }
  244.  
  245.     double IncompleteBeta(double x,double a,double b)
  246.     /*
  247.     The incomplete beta function Ix(a,b).
  248.     The Numerical Recipes routine assumes a>0. This work-around is taken from
  249.     Abramowitz and Stegun page 944, Eq. 26.5.16.
  250.     */
  251.     {
  252.         if(a>0.0)return betai(a,b,x);    /* Numerical Recipes in C */
  253.         else return IncompleteBeta(x,a+1.0,b)
  254.             +exp(gammln(a+b)-gammln(a+1.0)-gammln(b)+a*log(x)+b*log(1.0-x));
  255.     }
  256.     
  257.     double InverseIncompleteBeta(double p,double a,double b)
  258.     /*
  259.     This simple-minded bisection routine is slow, but its answer is accurate
  260.     to within ±1e-10. 
  261.     */
  262.     {
  263.         double f,low=0.0,high=1.0,mid;
  264.         int i;
  265.         
  266.         for(i=0;i<30;i++){
  267.             mid=(low+high)*0.5;
  268.             f=IncompleteBeta(mid,a,b);
  269.             if(f>p)high=mid;
  270.             else low=mid;
  271.         }
  272.         return (low+high)*0.5;
  273.     }
  274.     
  275. #endif
  276.